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Abstract: We propose a framework for solving the Einstein equation for static and Eu- 
clidean metrics. First, we address the issue of gauge-fixing by borrowing from the Ricci-flow 
literature the so-called DeTurck trick, which renders the Einstein equation strictly elliptic and 

00 generalizes the usual harmonic-coordinate gauge. We then study two algorithms, Ricci-flow 

and Newton's method, for solving the resulting Einstein-DeTurck equation. We illustrate 
the use of these methods by studying localized black holes and non-uniform black strings in 

qs, five-dimensional Kaluza-Klein theory, improving on previous calculations of their thermody- 

namic and geometric properties. We study spectra of various operators for these solutions, in 
particular finding the negative modes of the Lichnerowicz operator. We classify the localized 

^ solutions into two branches that meet at a minimum temperature. We find good evidence for 

a merger between the localized and non- uniform solutions. We also find a narrow window of 
localized solutions that possess negative modes yet have positive specific heat. 
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1. Introduction 

In four-dimensional general relativity, the space of physically relevant geometries, and the 
methods for finding them, are well understood. The stationary solutions are privileged, being 
possible end states of dynamical processes, and in vacuum the uniqueness theorems imply that 
the most general solution is the Kerr black hole, whose solution is known analytically. The 
problem of constructing dynamical spacetimes is in general a numerical one, as exact solutions 
are often unknown for situations of interest, in particular those with few symmetries such as 
generic gravitational collapse or black hole mergers. The tools for solving the dynamical 
Einstein equations are well understood and have enjoyed intense study both theoretically and 
computationally. 

The consideration of general relativity in more than four dimensions can be motivated 
in a variety of ways. Such studies date back almost a century to Kaluza and Klein, and 
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in a modern context are given weight by string theory and other phenomenological and 
formal considerations, such as Randall-Sundrum and large-extra-dimensions scenarios [1-4] 
and the AdS-CFT correspondence [5]. However the space of static and stationary solutions 
is considerably more complicated in these contexts, and many solutions known or believed to 
exist currently have no analytic solution. 

For asymptotically flat D-dimensional gravity in the absence of matter, with D > 4, static 
black holes are given uniquely by the Schwarzschild-Tangherlini black hole [6-8]. However, 
for stationary solutions in D = 5, in addition to the Myers-Perry solution which generalizes 
Kerr [9], the Emparan-Reall black ring is found [10]; these may even be combined into a 
'black Saturn' [11, 12]. Whilst these are understood analytically, in D > 5 the black ring 
solution is not known and it is thought that the structure of stationary solutions is yet more 
complicated than in D = 5 [13, 14], with little hope of admitting closed-form solutions. 

Now consider Kaluza-Klein theory, by which we mean gravity in five dimensions such 
that the asymptotic geometry is a product of four-dimensional Minkowski space and a circle. 
This is the canonical toy model for theories with compact dimensions. In this case the space 
of static vacuum solutions, which need not have an isometry in the compact direction, is not 
analytically understood, and the solutions that have been found numerically so far have a rich 
structure [15-17]. 1 A numerical approach exists in the literature to find cohomogenity-two 
static vacuum metrics, including black holes, employed first in the context of brane worlds 
[29] and then Kaluza-Klein theory to find non- uniform and localized black holes [30-32]. A 
similar approach was used independently in four-dimensional Einstein- Yang-Mills systems in 
the earlier work of [33,34] (including extension to four-dimensional stationary systems [35]) 
although the 'constraint equations' and their consistency was not discussed. The method 
has been used effectively in a variety of contexts [36-43], but is fundamentally limited in 
that it requires cohomogenity two, a particular conformal coordinate system, and a certain 
'wizardry' in actually making it work — for example, the 'trick' employed when treating a 
symmetry axis [15,29]. 

In Kaluza-Klein theory, the static vacuum solution that preserves four-dimensional Poincare 
invariance is trivial to write down analytically. However with additional compact dimensions 
even this is not true. For example, in string theory phenomenology the vacua of the the- 
ory that preserve supersymmetry typically involve products of four-dimensional Minkowski 
space with complicated manifolds such as Calabi-Yaus. In the supergravity approximation, 
and in the simplest cases with trivial fluxes and dilaton, the metric on the compact Calabi- 
Yau should be Ricci-flat, and such metrics are not known analytically. In general, finding 
four-dimensional Poincare invariant vacua can be thought of as the problem of solving the 

In the simplest examples of the AdS-CFT correspondence the situation is the same, with the 5-sphere 
playing the role of the Kaluza-Klein circle, and an analogous structure of solutions is thought to exist [18, 19]. 
The structure of static solutions here is particularly interesting, as the phase diagram of the black holes is 
dual to the phase diagram of the finite temperature gauge theory, and gravitational phase transitions can be 
translated into field theory ones [20-25]. Even black holes that are homogeneous on the 5-sphere may have 
a complicated form, such as the 'plasma-ball' black holes in simple confining models of AdS-CFT [26], which 
have interesting links to formalizing the membrane paradigm [27, 28] . 
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Euclidean Einstein equation on the compactification manifold with specific matter. The 
Euclidean problem of finding Calabi-Yau metrics that are Ricci-flat has received attention 
recently, with an approach based on Ricci-flow [44] and a spectral approach based on an alge- 
braic representation of the metric [45-48]. 2 However, these are fundamentally tied to working 
with Kahler geometries, and using the associated complex coordinates. 

A natural question is then whether there exist numerical methods to find general sta- 
tionary and static black holes in more than four dimensions, and to find solutions to the 
Euclidean Einstein equation such as those governing phenomenological vacuum compactifica- 
tions in string theory. In these higher dimensional settings the dynamical numerical methods 
that are so well explored in four dimensions simply generalize. Thus in principle a numerical 
algorithm exists: one chooses some initial data whose late-time evolution asymptotes to our 
stationary or static vacuum solution of interest. The first point to make is that this will never 
allow one to find dynamically unstable solutions, which may be of interest for theoretical 
reasons. However, the second more serious obstruction to this being a practical method is 
that dynamical evolutions, particularly involving horizons, are very complicated, partly due 
to the constraints that must be initially solved — already a hard elliptic problem in general — 
and then ensuring that these constraints are properly preserved under evolution. For these 
reasons we believe that such a dynamical approach is not a viable one in practice. 

The challenge then is to find a more direct but unified method to solve the Einstein 
equation in stationary, static, and Euclidean contexts. In this paper we provide a partial 
answer to this challenge. We present a general covariant framework in which to address the 
problem of finding static and Euclidean vacuum solutions in arbitrary dimensions, and with 
no particular symmetry properties (i.e. any cohomogenity) . The geometry may contain no 
horizon, or one horizon component, in which case the algorithm only yields the spacetime 
exterior to the horizon. We employ the fact that such static Lorentzian spacetimes, exterior 
to a horizon, may be analytically continued in time to smooth Euclidean geometries. The 
Lorentzian Killing vector field continues to a vector field which generates a U(l) isometry, 
and vanishes only at the point that continues to the Lorentzian horizon. Thus the problem 
of finding static vacuum solutions without horizon, or exterior to a single horizon, is included 
in the general Ricci-flat smooth Euclidean problem. The heart of the paper is to provide 
an algorithm to find smooth Ricci-flat Euclidean geometries in a fully covariant manner by 
directly solving the Ricci-flatness condition. 

We begin, in Section 2, by finding a general and covariant method for gauge-fixing that 
generalizes the harmonic-coordinate gauge condition. Specifically, we borrow from the study 
of Ricci-flow the so-called DeTurck trick [51], in which, using a fixed reference connection, the 
Ricci-flow equation is rendered strictly parabolic without changing its geometric content. In 
our context the DeTurck trick renders the Einstein equation into a strictly elliptic system of 
partial differential equations, which we refer to as the Einstein-DeTurck equation. In Section 
3, we then consider two types of algorithms for solving this elliptic system. 

2 The related case of Einstein metrics on del Pezzo surfaces has also been studied [49,50]. 
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The simplest algorithm is to exploit the ellipticity of the Einstein-DeTurck equation and 
solve it via a diffusive flow. The canonical flow is Ricci-DeTurck flow (diffeomorphic to Ricci- 
flow), and the idea is to provide an initial guess and then simulate the flow, which hopefully 
asymptotes to the desired Ricci-flat fixed point. However for vacuum black hole spacetimes 
the Ricci-flat solution may not be an attractive fixed point of Ricci-flow, as such black holes 
usually possess Euclidean negative modes [52] — in vacuum we expect that any black hole with 
a negative specific heat capacity is likely to have one or more negative modes [53-57]. Other 
interesting Euclidean geometries also possess negative modes, and show analogous behaviour 
to black holes under Ricci-flow [58]. We therefore present an algorithm in which a finite- 
dimensional family of initial guess geometries must be chosen, and the parameters of this 
family varied in order to reach the fixed point. Since the Ricci-DeTurck flow is diffeomorphic 
to Ricci-flow, which provides a flow not only in the space of metrics but also in the space 
of diffeomorphism classes (or 'geometries'), this method has the attractive feature that the 
algorithm is independent of the choice of reference connection. We denote this approach the 
'Ricci-flow method'. 

A more straightforward approach to solve an elliptic PDE is to iteratively linearize it, 
i.e. apply Newton's method. The point here is that there exist direct methods to solve 
the resulting linear PDE. For example, using real-space finite differencing this PDE appears 
as a sparse system of linear equations, for which fast solvers exist such as the biconjugate 
gradient method. Such solvers are not affected by negative modes (although they are by zero 
modes). Whilst more complicated to implement in practice, Newton's method avoids having 
to tune families of initial data — just one guess metric in the basin of attraction of the desired 
solution is sufficient. Unlike the Ricci-flow method, the trajectory of Newton's method in the 
space of geometries does depend on the reference connection, and hence so does the basin 
of attraction of a solution. In principle this makes the method somewhat less elegant than 
Ricci-flow, although in practice it appears to make little difference. 

In summary, finding Euclidean Ricci-flat geometries and static vacuum Lorentzian solu- 
tions, including black holes, can be phrased as the same problem and solved using the same 
unified, covariant elliptic approach. The problem of dealing with Euclidean negative modes 
in either context can be overcome either by using Ricci-flow and tuning initial data, or by 
using Newton's method. In the Lorentzian context, solutions may be found independently of 
their dynamical stability properties. Whilst these methods do not assume any isometries, if 
such isometries exist, it is natural to exploit them to effectively reduce the dimensionality of 
the problem. 

It is instructive to see the two methods described in Section 3 demonstrated in a non- 
trivial example. Thus in Section 4 we test both these methods by finding the various static 
black holes in five-dimensional Kaluza-Klein theory — the 'localized black holes' and 'non- 
uniform black strings'. We find that Newton's method gives excellent results in practice and 
is considerably easier to use than the Ricci-flow method with its tuning of initial data. We 
compute various properties of these black holes, such are their thermodynamic properties and 
horizon geometry. We are able to provide good evidence for a merger between the localized 
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and non-uniform solutions [15-17], computing the localized branch closer to the potential 
merger point than previously done in [59] even though we have used more modest resources. 

New results include finding that the localized solutions have a minimum temperature 
which divides the localized solutions into two branches which we term the small and large lo- 
calized black holes. 3 The Euclidean negative modes of the localized and non-uniform solutions 
are computed in the class of eigenfunctions which preserve the isometry of the background, 
and we find small localized solutions have a single negative mode, and the large localized and 
non-uniform strings computed have two. These negative modes are consistent with merger of 
the large localized and the non-uniform string branches, with one eigenvalue remaining finite 
and the other apparently diverging to minus infinity at the merger point. We also provide 
additional evidence for merger by considering low lying positive eigenvalues of various opera- 
tors. Perhaps the most intriguing result is the observation of a rather narrow window of small 
localized solutions that appear to have positive specific heat capacity, but possess a negative 
mode. This will require more rigorous numerical work to confirm. Assuming it holds true, 
this provides the first static counter-example to any claim that one might be able to reverse 
the link between local thermodynamic stability and the existence of negative modes discussed 
in [53-57], complementing the interesting recent stationary counter-example of [62]. 

2. The Einstein-DeTurck equation 

We wish to solve the vacuum Einstein equation, 



for a Euclidean-signature metric on a given D-dimensional manifold M, with appropriate 
boundary conditions if M has a boundary. For simplicity and concreteness, we are restricting 
our attention to the vacuum Einstein equation with zero cosmological constant. However, 
most of our considerations will be local, and including a non-zero cosmological constant and 
matter fields should not present any significant additional difficulties. 

As we will explain below, equation (2.1) is weakly elliptic; roughly speaking, it is elliptic 
for the physical degrees of freedom in the metric. Most methods for solving elliptic PDEs are 
relaxation methods, meaning that one starts with an initial guess for the unknown function, 
which presumably does not solve the PDE, and then iteratively deforms it to obtain functions 
that solve it to a better and better approximation. This is to be contrasted with methods for 
solving initial-value problems, associated with hyperbolic or parabolic PDEs, which usually 
proceed by building up the solution time-slice by time-slice. 

Before being able to apply any of the standard numerical methods for elliptic PDEs, one 
must deal with the issue of gauge invariance. The problem is that (2.1) is underdetermined. 

3 This nomenclature follows that of black holes in AdS [60] or in a cavity [61] which also possess a minimum 
temperature. However the thermodynamic behaviour of the localized black holes is rather different from these 
cases, with the large localized solutions having negative specific heat capacity unlike for AdS or a cavity. 
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Naively it appears to consist of D{D + l)/2 equations, per point, for D{D + l)/2 unknowns 
(the metric components). However, as a consequence of diffeomorphism invariance, the Ricci 
tensor for any metric satisfies the Bianchi identity, W^R^ — \d u R = 0, which has D compo- 
nents. Hence (2.1) actually represents only D{D — l)/2 non-trivial equations. Stated more 
precisely, the principal symbol of the Ricci tensor (considered as a nonlinear operator on 
the metric) annihilates pure gauge modes of the metric (corresponding to small diffeomor- 
phisms), while it is positive definite on physical modes. This is the technical meaning of 
"weakly elliptic". 

The same issue, of course, arises in the Lorentzian context, where the Einstein equation 
is only weakly hyperbolic, and our discussion will to a large extent parallel the standard 
treatment of the latter problem using harmonic coordinates (see for example section 10.2 
of [63]). As we will explain in subsection 2.1, the weak ellipticity makes the Einstein equation 
unsuitable for numerical solution without some form of gauge fixing. Of course, there are 
many different possible gauge-fixing schemes, many of them adapted to particular problems. 
Here we are looking for one that can be applied as generally as possible. A local analysis 
of the problem of gauge-fixing leads us to a generalized form of harmonic gauge that is 
somewhat more covariant than the textbook version. Then in subsection 2.2 we advocate an 
a posteriori method of implementing the gauge-fixing, in which the gauge-fixing condition is 
not imposed from the outset but rather solved for at the same time as the Einstein equation. 
Specifically, we propose solving not (2.1) itself but another PDE (see (2.11)) that we call 
the Einstein-DeTurck equation. The Einstein-DeTurck equation simultaneously encodes the 
Einstein equation and the gauge-fixing condition, and is strictly elliptic. In the next section, 
we will discuss algorithms for solving the Einstein-DeTurck equation. 

2.1 Gauge condition 

If we make a small perturbation bg^ = hu, v to the metric, the change in the Ricci tensor is, 
to first order, 

5R^ V = A R h^ u = A L h^ u + V^Uj,) , (2.2) 

where 

A L V = " 2 y2 V " R^u X h KX + R ifl K K )K , Vfl = V v h\ --d„h. (2.3) 

For the moment we are concerned with the local properties of Ajj, so we consider perturbations 
with wavelength much smaller than any curvature scale of the metric, and therefore take g^ u 
to be constant; this yields the principal symbol ^P rinci P al (the part of R^ u with only second 
derivatives of the metric) : 

A™ a V = -\^ 2 Ku ~ £ W + dxd^ . (2.4) 

Any such perturbation can be decomposed as h^ u = h^+duw^, where d v h u ^— ^d^h = 0; 
is the transverse (physical) part, and di^w v \ is the longitudinal (pure gauge) part. Specifically, 
the perturbation dr^w^ is the change in the metric under the small diffeomorphism generated 
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We will now briefly explain why the weak ellipticity of (2.1) is problematic for numerical 
purposes. For concreteness, suppose that we represent the metric by the values of its com- 
ponents on a set of lattice points in each coordinate patch (the argument does not change 
if we use some other representation). In general, to be representable the metric components 
must be smooth on the (coordinate) scale of the lattice spacing. However, this smoothness 
condition is not preserved by diffeomeorphisms. Thus, even if our lattice is sufficiently fine to 
represent some particular metric that solves (2.1), other metrics in its difieomorphism class 
will not be representable. (In some sense, "most" such metrics will not be.) Therefore we need 
to control which representative of the difieomorphism class we aim for. What will go wrong 
in practice if we do not? During the relaxation process, since short-wavelength pure-gauge 
modes are not damped out, small errors will accumulate, leading to a loss of control. 

There are of course many different ways to fix a gauge in general relativity. One can fix 
certain metric components; for example, in the presence of spherical symmetry one can set the 
coordinate r to equal the proper radius of the sphere. Alternatively, one can take advantage 
of certain special properties of the metric. If one is solving for a Kahler metric, then one 
can choose to employ complex coordinates; the only allowed gauge transformations are then 
holomorphic diffeomorphisms, which have no local degrees of freedom. (This is actually a 
special case of harmonic coordinates, which we will discuss below.) Our aim in this paper, 
however, is to find an all-purpose method for solving (2.1), one that requires as little special 
knowledge of the metric one is trying to find as possible. 

Let us return to our local analysis. We wish to forbid the pure-gauge fluctuations; in 
other words we wish to require d v h v ^ — \d^h = 0. We need to de-linearize this equation, i.e. 
to find a one-form depending on the metric (and of course not diffeomorphism-invariant) 
whose variation 5^ about a constant equals d„h v ^ — ^d^h. Fixing £^ will then project 
out the pure-gauge modes. The simplest such one-form is 



where are the coordinates considered as scalar functions and V 2 is the scalar Laplacian. 
The gauge choice £ M = is the harmonic gauge, commonly employed in Lorentzian numerical 
relativity. 

Although it does the job, the harmonic gauge is a bit inelegant insofar as it refers to 
a particular coordinate system. In particular, transforms in a complicated way on chart 
overlaps, and a metric that is in harmonic gauge in one patch will generally not be in that 
gauge in a neighboring one. This problem is easily fixed by replacing the coordinate derivatives 




(2.5) 



Here are three expressions for the dual vector field 




(2.6) 
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in (2.5) with any fixed covariant derivative. We thus choose an arbitrary background metric 
g^y on M, and define 

= a"" " 2 ^a) , ^ = - f^) = - (f) V„ ( (|) 9 A , 

where f 1 ^ is the Levi-Civita connection for g^ and is the associated covariant derivative. 4 
In the theory of harmonic maps between Riemannian manifolds, ^ is referred to as the tension 
field of the identity map from (M,g^ u ) to (M,g^ v ), and the gauge-fixing condition 

e M = (2.8) 

is equivalent to the requirement that the identity map be harmonic. 

Ideally a gauge-fixing condition should intersect each gauge orbit exactly once: given a 
metric Qav, there should be a diffeomorphism (unique up to isometries) that maps it to a 
metric g^ v satisfying (2.8). Such a diffeomorphism can be composed with the above identity 
map from {M,g^, v ) to (M^g^) to obtain a harmonic map from (M, g^) to (M, g^v). In 
physicists' language, this harmonic map is a solution to the classical equation of motion for 
the sigma model from (M,g^ u ) into (M,g^ u ): 

S[y}= [ g^g^g^d^d^, (2.9) 

JM 

where y a are the coordinates on the target space (M, g^ u ). At the local level, it is easy to see 
the existence and uniqueness of solutions to (2.8). Under a small diffeomorphism, generated 
by a vector field w^, the change in is <5£^ = —Ayw^, where 

A v w» = ~V 2 ^ ~ \RV + (Tl - fl)V x w» . (2.10) 

Ay is a strictly elliptic operator, so locally £^ + 5^ = will have a unique solution for w^. 
The global existence and uniqueness problems are more difficult, of course. On a compact 
manifold the operator defined by (2.10) will have at most a finite-dimensional kernel, so the 
space of harmonic maps is at most finite-dimensional (see also [66]). As far as existence is 
concerned, there is no general proof, but there are many partial results (including for the case 
of manifolds with boundary); for example, Eells and Sampson [67] proved existence under the 
assumption of non-positive sectional curvatures for g^. As far as we are aware, there are no 
known counterexamples. 5 

4 Note that this is different from the gauge choice referred to as "generalized harmonic coordinates" in the 
numerical relativity literature [64,65], which involves setting V 2 a; M = fP , where is a fixed vector field. Of 
course one could combine the two generalizations by setting £ M = ff M . 

5 See [68] for a summary of the literature on harmonic maps through 1991. 
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2.2 A posteriori gauge fixing 

We now turn to the problem of implementing our gauge-fixing condition = 0. This is 
a differential equation for the metric, in contrast to algebraic gauge-fixing conditions such 
as fixing the values of certain metric components. In the latter case, when implementing 
a relaxation method, one may simply impose the gauge-fixing condition from the outset. 
Assuming the gauge fixing is complete, the Einstein equation will be a strictly elliptic operator 
on the space of metrics obeying the gauge-fixing condition. In our case, however, the space 
of gauge-fixed metrics is not so easy to describe explicitly. 

Instead of imposing the gauge-fixing condition a priori, we therefore propose to solve 
£^ = simultaneously with R^ u = 0. This can be done quite elegantly by combining the two 
equations as follows: 

<, = ^-V ( ^) =0, (2.11) 

which we will refer to as the Einstein- DeTurck equation. The tensor Rjf u appears in work of 
DeTurck [51]; he gave a simple proof of short-time existence for the Ricci-flow, dg^ u /d\ = 
—2R t j iU , by pointing out that it is diffeomorphically equivalent to the strictly parabolic PDE 
dg^y/dX = —2RH. Roughly speaking, the term V( M £j,) supplies the longitudinal components 
that are missing in the Ricci tensor. The linearization of R^ v is 

5R* V = A H V = A x V " \ £ i V + v (m*M > «" = " ^ uX W x • (2-12) 

The second and third terms on the right-hand side do not contain second derivatives of h^u, 
so the principal symbol of R^ is simply — ^ V 2 /^, and (2.11) is strictly elliptic. In effect, 
rather than projecting out the pure-gauge modes (which would be the result of working in 
the space of gauge-fixed metrics), we have given them a kinetic term, putting them on the 
same footing with the physical transverse modes. 

We now turn to the question, does a solution to (2.11) necessarily solve (2.1)? We begin 
with a local analysis. Equation (2.11) implies 

V^JL - \d v R H = ~\ (V 2 ^ + RS l Q = , (2.13) 

which locally has only the trivial solution £^ = 0. Thus Rjf u = 0, which is D(D + 1) /2 
equations for the same number of metric components, is locally equivalent to the D{D — 1) /2 
equations of R^ u = and the D equations of £ M = 0. 

As usual, the global question is more difficult. However, Perelman has shown that, on 
a compact manifold without boundary, (2.11) implies (2.1) [69]. (In Perelman's theorem, 
is a general one- form; it does not refer to the particular form (2.7).) It would be interesting 
to generalize his theorem to manifolds with boundary. (One would have to impose some 
boundary condition on such as the vanishing of its normal component. Otherwise one can 
easily construct a counterexample by taking any compact region of any of the known non- 
compact steady Ricci solitons.) The theorem also covers the case with negative cosmological 
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constant: R^ u — V = ng^, with k a negative constant, implies R^ u = Kg^. On the 
other hand, for positive k there are known counterexamples (called shrinking Ricci solitons), 
for example on the first and second del Pezzo surfaces [70-73]. However, the existence of 
such metrics should not present a major problem as long as they are isolated, since it is 
easy enough to check whether one's algorithm is converging to a true solution to the Einstein 
equation or merely to a soliton, and in the latter case to restart it with a different initial 
guess. (Similar comments presumably apply in the presence of matter fields, but we have not 
studied this case in detail.) 

Let us now discuss what boundary conditions are appropriate for equation (2.11). It is 
well known that the Einstein equation, being effectively D{D — l)/2 equations for the same 
number of variables, requires the same number of boundary conditions. For example, one 
can fix the induced metric on the boundary or its extrinsic curvature, which are respectively 
Dirichlet and Neuman boundary conditions for the tangential components of the metric. For 
the Einstein-DeTurck equation, we have an additional D equations and therefore need an 
additional D boundary conditions. Since we want to obtain a solution with £^ = 0, and since 
(2.13) is strictly elliptic, it seems clear that we should require ^{dM = 0. If we denote the 
normal and tangential components of £ M by £ n and respectively, then it is easy to show that 
£, n \dM = fixes the normal derivative of g nn , while C\dM = fixes the normal derivative of 
g n i, in terms of their values, and the values and derivatives of the other metric components. 
(This is the same issue as imposing boundary conditions for the parabolic PDE mentioned 
above, dg^/dX = -2Rf LU . In [52] it was argued that ^\om = are appropriate boundary 
conditions for that PDE, and uniqueness has been proved in that context [58].) 



3. Algorithms for solving the Einstein-DeTurck equation 

The Einstein-DeTurck equation (2.11) is a strictly elliptic nonlinear PDE, and there exist 
many standard methods for solving such equations. As mentioned above, most of these are 
relaxation methods, meaning that one starts with a trial function that does not solve the 
equation, and then iteratively deforms it to obtain functions that solve the equation to a 
better and better approximation. There are two basic strategies for how to relax the metric. 
If we iteratively replace g^ u with g^ + h^, then one can either deform in the direction of 

_r>H 

hfjtv = -tR^u , (3.1) 

or solve the linearized equation, 

V = -A^i?^. (3.2) 

The first is closely related to Ricci-flow, and the second is a Newton method. In this section 
we will discuss these two approaches in turn, trying to understand their general properties 
(rather than the specifics of how to implement them) . Our discussion should not be considered 
as exhaustive, and these strategies can be generalized, refined, and combined in various ways. 
We will for the most part work in the continuum limit, although it should be understood 
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that in practice the metric and operators will be discretized in some form. The choice of 
discretization scheme — finite differencing, finite element, spectral, etc. — is a separate issue 
that we will not address in this section. (For the application discussed in Section 4, we 
have chosen finite differencing, because it is the easiest to implement. One intriguing other 
possibility would be to implement a finite-element method on the basis of the Regge calculus.) 

Before proceeding, it will be useful to understand the spectrum of the operator Ajj — the 
linearization of R^ v — about a metric satisfying R^ v = and £ M = 0. 

3.1 Spectrum of Ah 

With = and £ M = 0, the operator Ah simplifies as follows; we will also need Ay (see 

(2.2) , (2.3), (2.10), and (2.12)): 

A H V = + v (mM . *" = - S Z" = + \ d ^ h + ( r ^A " Kx) h " X ■ (3-3) 

A v w* = -\v 2 w» + (T£, - f£, )VV . (3.4) 

Ar is a self-adjoint operator on the space Ti of metric fluctuations, with respect to the 
inner product j M g 1 l 2 (h tlv h IJ,u — \h 2 ), and therefore has a real spectrum. Let Tii C Ti be the 
space of pure-gauge fluctuations, i.e. those of the form h^ u = V^w^ for some vector field w^. 
Ti\ is annihilated by Ar, since the Ricci tensor remains zero under a gauge transformation. 
The quotient space Ti/Tii is the space of physical metric fluctuations. On this space Ar 
may have a finite number of negative or zero eigenvalues. For example, if there is a moduli 
space of solutions, then each modulus will give rise to a zero mode (although the converse 
does not hold: not every zero mode can be "exponentiated" to produce a continuous family 
of solutions). The Euclidean Schwarzschild metric has a single negative mode, concentrated 
near the horizon [74], as does the "small black hole" solution for gravity in a cavity [57]. 

We now turn to the spectrum of Ah- Note that this operator is not self-adjoint, and 
may have complex eigenvalues. It maps the pure-gauge fluctuation V^w^ to V ' ^Ayw u y It 
therefore has the same spectrum as Ay acting on the space of vector fields (modulo Killing 
fields), which may include a finite number of eigenvalues with negative or zero real part. From 

(3.3) we see that the action of Ah equals Ar when acting on the quotient H/Hi, so they 
have the same spectrum. 

In the presence of a boundary, there is one technical subtlety we should address. For 
concreteness we assume that we are imposing Dirichlet boundary conditions on the metric. 
This puts a boundary condition hij\gM = on the metric fluctuations in Ti. Furthermore, 
the vector fields w^ 1 generating elements of Hi should vanish on the boundary. Now, when we 
study the operator Ah, we impose an extra boundary condition, i^|aM = (the linearization 
of £^ = 0). Let Ti' C TC be the space of fluctuations obeying this boundary condition, and 
Ti\ = Tii H Ti' be the space of such pure-gauge fluctuations. If, as we have been assuming, for 
each metric there is a gauge transformation taking it into the gauge £^ = 0, then Ti' /Ti', = 

n/Hi. 
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All in all, the spectrum of Ah is the union of the spectrum of Ar (on TL/TLi) and the 
spectrum of Ay (on TL\). 

3.2 Ricci-flow method 

We consider equation (3.1). Here we are following the strategy of solving an elliptic PDE 
by simulating the associated parabolic one, for example, solving the Poisson equation by 
simulating diffusion. In the limit that we take the parameter e very small, we are simulating 
Ricci-DeTurck flow [51], 

^ = -2<. (3.5) 

The sign is chosen to make this a parabolic flow. Since locally the operator An is — ^V 2 , 
locally we are indeed simulating diffusion, and the same issues arise as in that case. For 
example, in order to avoid short-wavelength instabilities e must be chosen positive but not 
too large; the higher the spatial resolution the smaller e must be. 6 The method can be sped 
up by various strategies, such as Gauss-Seidel updating, successive over-relaxation, and the 
multi-grid method. 

An important property of the flow (3.5) is that it is equivalent, modulo a time-dependent 
diffeomorphism, to the pure Ricci-flow dg^/dX = —2R^ U . Therefore the diffeomorphism 
class of the metric evolves independently of the choice of background metric g^y. 

Numerical simulation of Ricci-flow has been performed previously in various contexts 
[49,52,58,73,75-78]. In [49] it was used as a method for solving the Einstein equation (with 
positive cosmological constant), on the third del Pezzo surface. There the flow was performed 
in the space of Kahler metrics, and a theorem of Tian-Zhu guaranteed that any initial Kahler 
metric (in the correct class) would flow to the Kahler-Einstein metric [79]. 

What can we say in the absence of such a theorem, for example in the case of black hole 
metrics? Supposing that a solution to R^ v = (and £ M = 0) exists on the given manifold with 
the given boundary conditions, can we expect Ricci-flow starting from a generic initial metric 
to converge to it? The necessary and sufficient condition for generic convergence is that the 
fixed point be attractive. The stability of a fixed point under the flow dg^y/dX = —2R^ V is 
controlled by the spectrum of the operator Ah, which we studied in the previous subsection. 
We found that the spectrum of Ah is the union of the spectrum of Ay (acting on the space 
of vector fields) and the spectrum of Ar (acting on the space of physical metric fluctuations). 
The former depends on the choice of background metric g^ u as well as g^y, while the latter 
depends only on g^ v . Either operator may have a finite number of negative eigenvalues (or, 
in the case of Ay, eigenvalues with negative real part). Such an eigenvalue for Ay would 
correspond to an instability in the gauge-fixing condition, and could presumably be cured by 
making a more judicious choice of g^ v . (For example, if g^ happens to equal the fixed-point 
metric g^y, then Ay = — ^V 2 , which has a positive spectrum.) In our experiments we have 
not so far encountered such an instability. On the other hand, negative eigenvalues of Al are 



6 This limitation can be circumvented by implementing an implicit updating scheme. However, this is 
probably not worth the trouble unless one is interested in accurately simulating the Ricci-flow per se. 



- 12 - 



an invariant property of the metric, which occur for many types of Euclidean black holes that 
one might wish to find numerically, including the ones we will study in Section 4. Naively, one 
might discard Ricci-flow as an unuseful algorithm for finding such solutions. However, we will 
now argue that, given enough information about the nature of the flows for different initial 
metrics, this problem can be overcome by supplementing the Ricci-flow with a root-finding 
algorithm on the space of initial data. 

Let us first set up some useful notation. Let M. be the space of diffeomorphism classes 
of metrics on M satisfying the given boundary conditions. (By a slight abuse of terminology, 
we will refer to these diffeomorphism classes as metrics.) Let g^ be the Ricci-flat metric we 
are searching for. The tangent space to M. at g( R ) is (isomorphic to) H/Hi- For simplicity we 
assume that there are no zero modes of A^ about g( R \ and in particular that there is not a 
moduli space of solutions. Let E be the set of points in M. through which a flow passes that 
asymptotes to g^ in the future (A — ► oo) , and T the set of points through which a flow passes 
that asymptotes to g( R ) in the past (A — > — oo). The tangent space to E at g^ is spanned by 
the positive modes of A^, and the tangent space to T by its negative modes. Hence, if N is 
the number of negative modes, then T is iV-dimensional and E is iV-codimensional. Starting 
from a point close to but not on E, the flow will pass close to g^ before being diverted by the 
negative modes and asymptotically approaching the surface J 7 . (The flow may alternatively 
hit the boundary of M at finite A due to the formation of a metric singularity, but it will be 
exponentially close to T when this happens.) 

The easiest case is N = 1. Then E is a hypersurface in M, and (at least locally) divides it 
into two regions, while T consists of two flows emanating in opposite directions from g^ R \ A 
concrete example is provided by the "small" Euclidean Schwarzschild black hole in a spherical 
cavity. The Ricci- flows of this system, for D = 4, were studied in the paper [52]. Let us briefly 
review that work. The small black hole is a Ricci-flat metric on the manifold B 2 x 5 2 , with a 
U(l) x 50(3) isometry group. Its negative mode shares those isometries and is localized near 
the horizon (the fixed locus of the U(l) action). In [52] the two flows emanating from the 
small black hole were found numerically. For one of them, the horizon shrinks to zero size in 
finite flow time, creating a singularity, while in the other the horizon expands monotonically 
and the metric asymptotically approaches the "large" black hole solution, which is of the 
order of the size of the cavity. The flow starting from any metric sufficiently close to, but not 
on, E will exhibit one of these two behaviors. Thus it is straightforward to determine — even 
automatically — on which side of the E surface a given metric lies. If the small black hole 
solution were not known, we could find it starting with any one-parameter family of initial 
metrics that intersects E by using a bracketing method to find that intersection point. 

For small localized black holes in Kaluza-Klein theory, as we will discuss later in detail, we 
expect a similar situation. There exists one negative mode, which either shrinks the horizon 
sphere (forming a singularity) or expands it (also forming a singularity since the analogous 
large black hole cannot be reached without a topology change). Thus one does not need to 
know the explicit form of the negative modes of A^. One can simply deduce which side of E 
a flow lies on, and use a bracketing procedure. 
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As we see later the non-uniform strings and some localized solutions have two negative 
modes. Then the situation is more interesting as the surface £ is codimension-2. Whilst here 
we have not attempted to solve this shooting problem, presumably it is still possible to do so, 
although it would be more complicated as one cannot use simple physical arguments to deduce 
what flow one is tending towards on the 2-dimensional surface T. One would first tune the 
initial data to remove the dominant negative mode (that with highest absolute eigenvalue), 
and then tune again to remove the second negative mode. 

There are (at least) three potential difficulties to be addressed in the discussion above. 
Firstly, the picture presented presupposes we can construct a family of initial data X that 
intersects S. Whilst an N dimensional surface may generically intersect a codimension N 
surface, it may also not. Thus one cannot naively take an N dimensional family of initial 
geometries. One must use some physical intuition to determine how these must vary in order 
for the family to contain a point in S. As we will see later in the Kaluza-Klein example this 
is certainly possible to do for one negative mode, but we have not tried for more negative 
modes, and presumably it becomes harder to find such data. Secondly, having found data 
I that intersects S, one must consider points on X that are sufficiently close that they are 
initially attracted to the fixed point of interest. Points far from £ may flow directly to other 
fixed points or to the boundary of A4 (by formation of singularities). In any algorithm that 
requires an initial guess there is always a basin of attraction around a fixed point, and the 
difficulty here is that one must lie in the basin of attraction for a multi-parameter set of 
initial data, rather than just finding one point in that basin. The third difficulty is that for 
a complicated black hole we may not know a priori how many negative modes exist. This 
is a more severe problem. One approach would be to start with a one-dimensional family 
of initial data, and if that doesn't work move to two and then three. However, given the 
above-mentioned problem that without knowing how to intersect £ one might simply miss it, 
it seems unlikely that one could use this method in practice without some physical intuition 
on the number or nature of the negative modes. 

3.3 Newton's method 

We now turn to (3.2), in which we iteratively solve the linearized equation. Actually, it is 
useful to refine (3.2) a bit to allow for a variable step size: 

V = -eA^fl* , (3.6) 

where e is a parameter (generally between and 1) that may be chosen independently at 
each step. This is useful because otherwise — particularly if the inital metric is not close to 
a solution — some steps may be extremely large, and in fact may lead outside the space of 
(positive-definite) metrics. How to choose e optimally is a standard problem in the imple- 
mentation of Newton's method; in fact for the Kaluza-Klein examples we exhibit later we 
were always able to take e = 1. 

The utility of the method lies in the fact that there exist direct methods, such as bicon- 
jugate gradient, to solve the linear equation A^/i^ = —eR^. In a real-space discretization, 
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the operator Ah is sparse, and these methods are extremely fast. They work well as long 
as Ah is well conditioned (does not have eigenvalues that are very small in absolute value). 
In contrast to the Ricci-flow method, they are not in the least affected by the presence of 
negative eigenvalues. Hence any fixed point will possess a genuine basin of attraction in the 
space of initial metrics — no root-finding algorithm is required. 

Another difference with the Ricci-flow method is that, in (3.6), depends on the choice 
of background metric j™, in a way that cannot be removed by a diffeomorphism. In this sense 
the method is less "inherently geometrical" than Ricci-flow. The latter describes a continuous 
curve in A4 (the space of diffeomorphism classes on M) depending only on the initial metric, 
with an evolution that is local on M. Newton's method, on the other hand, jumps around 
in the space of metrics in a way that is neither diffeomorphism-invariant nor local. Thus, 
unlike Ricci-flow, it neither requires nor provides information on the global structure of A4. 
It simply homes in on the nearest fixed point. In practice we found it to be a more robust 
algorithm than Ricci-flow in the case where negative modes of occur for the solution of 
interest. 

Ricci-flow and Newton's methods also differ in their rate of convergence. Whilst the 
approach to a fixed point under Ricci-flow is governed by the spectrum of the operator A^, 
with short wavelengths converging in flow time quicker than long, for Newton's method the 
convergence in iteration time is uniform in all wavelengths. This is a potential advantage of 
the method when considering non-compact manifolds, where one might truncate the manifold 
for numerical purposes, but then removing this regulating truncation A l will have modes with 
increasingly small eigenvalue which would converge very slowly in Ricci-flow time. 7 

4. Application to black holes in Kaluza-Klein theory 

We will now use the example of black holes in 5-dimensional Kaluza-Klein theory in order to 
illustrate the two algorithms we have presented here to solve the Einstein-DeTurck equation. 
To this end we begin with a brief review of the various black hole types we consider, before 
moving on to show how to apply the methods and present results. 

4.1 Review of black holes in Kaluza-Klein theory 

We will be concerned only with single component horizons. The geometries of interest will 
asymptote to the flat Euclidean metric with compactified Euclidean time r, with proper 
length j3, and Kaluza-Klein circle direction with length L. Thus the geometry asymptotes to 
the product x R 3 x S\ where Sp is the Euclidean time circle and is the Kaluza-Klein 
compact dimension, with metric, 

ds 2 = dr 2 + dy 2 + y 2 dn 2 + dz 2 (4.1) 

7 Note that we are referring to Ricci-flow time and the iteration time for Newton's method, rather than the 
actual computer time. For example, in practice the time taken to solve the linear system required for each 
iteration of Newton's method may depend on the condition number of the linear system which will become 
worse if very small eigenvalue modes occur. 
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using angular coordinates r G [0, (3) and z £ [0, L), and we have used spherical coordinates on 
M 3 with dS7 2 the line element on a unit S* 2 . Since we are interested in static black holes then 
d/dr will be a Killing vector field and generate a U(l) isometry, and will vanish smoothly at 
the horizon. 

The known black holes preserve the 50(3) axisymmetry manifest in the asymptotic 
geometry above, and thus to a lower dimensional observer they appear spherically symmetric. 
Hence they can be written in the form, 



where T, A, B, F, S are functions of two coordinates x, y. The solutions fall into 3 classes; 
Localized black holes 

This one parameter family of solutions have a horizon topology S" 3 . For small horizon area 
the geometry near the horizon appears similar to the 4-dimensional Schwarzschild solution. 
As the area is increased the horizon geometry deforms as it 'feels' its images in the extra 
dimension. Little is known about this branch of solutions in the highly deformed region. 
Previous numerical constructions for D = 5,6 have been given by [31,32,59], and analytic 
approximations have been used to study the branch in the small area limit as an expansion in 
the horizon radius over L [80-85]. We find later that the small localized black holes possess 
one negative mode of as for 4-dimensional Schwarzschild. As one moves along the branch 
of solutions, increasing their size, their temperature decreases. We term these solutions the 
small localized solutions. A minimum temperature is then reached. At this point the tangent 
to the solutions gives a zero mode of A^. Proceeding further the solutions become hotter 
again, and the zero mode continues to a second negative mode. We term these solutions the 
large localized solutions. 

Uniform black strings 

These are a one parameter family of solutions with horizon topology S 2 x S 1 which we may 
write explicitly as the product of 4-dimensional Schwarzschild with a circle, length L; 



where ro provides the parameter along the family. Due to the direct product structure of the 
metric these solutions inherit the negative mode of the 4-dimensional Schwarzschild solution. 
Using the above coordinates the Schwarzschild negative mode is a transverse traceless mode of 
A^ preserving the SO(3) isometries of the form, A^\r), with eigenvalue — oJ 2 r 2 ~ 0.38 [74]. 
Then the negative mode inherited by the uniform string takes the form h^ u = Afy\r), 
hfiz = hzz = 0, and has the same eigenvalue, so that it preserves all the isometries of the 
uniform string, including the translational invariance generated by d/dz. 

Gregory and Laflamme demonstrated that for a critical ro which we call r c , this solution 
has a zero mode of the Lichnerowicz operator [86] . Reall showed that this is simply related to 



ds 2 = Tdr 2 + Adx 2 + Bdy 2 + 2Fdxdy + Sdfl 2 



(4.2) 



ds 2 = (1 - r /r) dr 2 + (1 - ro/r)" 1 dr 2 + r 2 dtt 2 + dz 2 



(4.3) 
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the usual 4-dimensional Schwarzschild negative mode [56] with the Gregory-Laflamme zero 
mode of the uniform string with radius r c being given by A^\r)e luJZ . For tq < r c we expect 
this zero mode to flow to a negative mode with similar harmonic z dependence, and this 
is indeed what we see numerically - in Lorentzian signature the zero mode extends to the 
dynamic Gregory-Laflamme instability. Conversely for ro > r c the mode flows to a positive 
mode. Hence the uniform strings fall into those with ro > r c and one translationally invariant 
negative mode, and those with ro < r c for ro near r c which have in addition a negative mode 
with harmonic z dependence. In fact, as one decreases ro further one obtains additional 
higher harmonic negative modes, associated with higher harmonic zero modes. 

Non-uniform black strings 

The zero mode of a uniform string with ro = r c provides a linear Ricci-flat deformation of the 
background solution. This linear deformation can be lifted to a full non-linear deformation, 
as illustrated by Gubser's order by order analysis [87]. The new branch of solutions, which 
is inhomogenous in the circle direction, are hence termed non-uniform strings. Numerical 
constructions of this branch can be found in [30,39,40]. We find later that the non-uniform 
strings emerging from the uniform branch have two negative modes. One is the continuation 
of the homogeneous uniform string negative mode. The second is the continuation of the zero 
mode at the branching point into a negative mode, agreeing with the Morse theory picture 
of Kol [15,88]. 

The uniform strings are the only solution known analytically, and therefore the task is 
to employ our numerical method to construct the localized and non-uniform string solutions 
about which relatively little is known. A natural conjecture is that the localized black holes 
and non- uniform string branches meet at a 'merger point' [88-91]. In particular Kol has 
made a specific proposal for the singular solution that enables the change in topology [88]. 
Studies of the branches from both side indicate consistency with such a 'merger' [59,92], and 
local studies of the increasingly curved region of the geometry that is supposed to become 
singular also indicates consistency [39,93]. However there are, as yet, no conclusive studies, 
largely because the localized branch of solutions has proven more difficult to construct than 
the non-uniform solutions which have been recently constructed to high precision [39] . 

4.2 Characterizing the Kaluza-Klein black holes 

We consider our static vacuum black hole spacetimes as smooth Riemannian manifolds with 
a U(l) Killing vector that vanishes at the horizon. Then we are naturally led to work with 
the canonical ensemble for the black hole, with the black hole Hawking temperature T = 1/(3 
being the control parameter, or classically, the asymptotic size (5 of the circle generated by 
the Killing vector. In addition to this control parameter we have the asymptotic Kaluza-Klein 
circle size L. Due to the classical scale invariance of the vacuum Einstein equations we may, 
without loss of generality, choose to fix L and work in units where L = 1 . Solutions with any 
other value of L may then be trivially obtained by the appropriate scaling, 
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Since the static black hole geometries are smooth with only an asymptotic boundary, then 
fixing L we have that j3 is the only control parameter available. However previous numerical 
results [59] indicate (although the effect appears subtle so they are inconclusive) that the 
localized black holes have a minimum temperature as one moves along the one parameter 
branch of solutions. This is indeed confirmed by our numerical results here. For generic 
static black hole solutions one might expect that branches of solutions have maxima and 
minima of the temperature. This raises two important questions. Firstly how does one find 
the solution of interest given that two or more solutions may have the same topology and 
temperature, ie. there is no black hole uniqueness? Secondly, suppose in a given topology 
there are no solutions below a certain temperature - as we will find later for the localized 
black holes. Then how do we see that no solution exists? 

The answer to these questions depends on the method being used. Let us first discuss 
Newton's method. For a given temperature multiple solutions may exist. However, the basis 
of attraction of the two solutions will be different. Hence one must find an appropriate initial 
guess geometry that selects the solution of interest by being initially sufficiently close to it. 
We will see this in practice later for the localized black holes. For Ricci-flow, black hole fixed 
points may not be attractors due to negative modes. However qualitatively the same answer 
is true, so that after appropriately tuning initial data to arrive at a fixed point, the family of 
initial data we tune will determine the fixed point reached. 

If no solution exists at a given temperature then Newton's method will find no solution. 
However, as we have discussed above, Newton's method does not give a flow in the space 
of geometries (diffeomorphism equivalence classes) so if one does not find a solution one 
cannot simply conclude that one does not exist. Since the Ricci-flow is a flow in the space 
of geometries, one might hope that its global behaviour does give an indication of whether 
a solution exists or not. Unfortunately this is not so clear cut as the tuning of the family of 
initial data requires that that family intersect the finite codimension hypersurface £ which 
is closed under Ricci-flow and attracted to the fixed point. Failure to find a solution can 
simply indicate that as the control parameter changes, the family of initial data chosen - now 
a function of the tuning parameters and this control parameter - fails to intersect S. 

Thus we control which black hole we wish to locate by choosing the appropriate topology 
(say to differentiate strings from localized black holes) and then specify the inverse temper- 
ature (3, and having done that choose an appropriate initial guess. Having found solutions 
there are various quantities of interest that we will display later. 

Firstly at large radius, S ~ y 2 — > oo, the time and Kaluza-Klein circle sizes become 
homogeneous in x, and these homogeneous components give rise to 2 asymptotic charges, 
the mass (which is conserved under time evolution in the Lorentzian setting) and relative 
tension (which is not conserved). Following [94-97] this leading homogeneous correction to 
the asymptotic metric takes the form, 




(4.4) 
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choosing the radial coordinate as the sphere radius. Then, choosing units where G/v = 1 the 
mass M and tension r are defined as, 



M 




t L = - a — b 
2 



(4.5) 



and give rise to a first law and Smarr relation, 



dM = TdS + rdL 



2M = 3 TS + tL 



(4.6) 



where S, the entropy, is given in our units by S = Ah or i z /A for horizon area Ah or iz- As usual 
the free energy is given by F = M — TS. 

Secondly the geometry of the solutions, particularly near the horizon is of interest. For 
the localized solutions we define R eq to be the round equatorial 2-sphere radius of the horizon, 
and Lp i ar to be the geodesic distance from the equator to the pole. The axis of symmetry 
is exposed in the localized solutions, and we define 2L ax i s to be the proper length of the 
geodesic between the poles of the horizon that follows this axis about the Kaluza-Klein circle. 
We define the eccentricity of the horizon, e, to be given in terms of the ratio of the area 
of the round equatorial horizon 2-sphere, A eq and the geodesic 2-sphere in the horizon that 
contains both poles A po i, as e = A po i/A eq — 1. We denote R T to be the proper length divided 
by (27r) - the 'radius' - of the time circle at the point on the symmetry axis equidistant from 
the horizon poles. Finally we may regard the time circle as being fibered over the axis of 
symmetry to give s 2-sphere. We call this the time 2-sphere, and call its area A^ me . For the 
non-uniform solutions we denote Rmaxi Rmin the maximum and minimum round 2-sphere 
radii within the horizon - for the uniform strings we take these both equal to the horizon 
radius. We define 2Lf lor i z to be the proper length of a closed geodesic within the horizon that 
wraps the Kaluza-Klein circle once. 

We will also consider the geometry of the horizon directly by embedding it into Euclidean 
space later. However the geometry near the horizon is also of interest, and one probe of this 
that we will consider is the eigenvalue of the negative modes of that the metric possesses. 

4.3 Details of numerical construction 

We now discuss some pertinent details that underlie the numerical construction, and are 
important in order to correctly interpret the results. 

We use second order accurate real space finite differencing to approximate the Einstein- 
DeTurck equations. In order to implement the Kaluza-Klein asymptotics we introduce a 
boundary at large but finite 2-sphere radius and fix the induced metric on the boundary to be 
that of the round sphere, radius R » L, with thermal circle length (5 and Kaluza-Klein circle 
with length L. Thus we are making two approximations, firstly the finite differencing and 
secondly the finite asymptotic boundary. In appendix A we discuss convergence testing and 
dependence of our results on the physical boundary position, justifying the approximations 
used here. In addition we demonstrate that whilst we are solving the Einstein-DeTurck 
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Figure 1: The geometry of the two overlapping charts for the localized black holes. The points 
drawn are for the lowest resolution considered here, 20 x 60 with R = 3L where we have chosen L = 1. 
Boundary points are shown enlarged as are points that must be interpolated from the other chart. 

equations, the vector field £ is vanishing in the solution which is, therefore, Ricci flat as 
required. 

For the non-uniform strings we may choose one coordinate chart with coordinates as 
above (4.2), where we take x to be the circle direction and y the radial direction. Then 
y = and R give the horizon and truncated asymptotic boundary, and x = and L have 
mirror boundary conditions to implement the periodic Kaluza-Klein circle. It is an important 
point that the generalized harmonic coordinates we use have sufficient freedom to choose the 
boundaries of the chart to be at these locations. For the localized black holes we use two 
overlapping coordinate charts, one adapted to the asymptotic region, given by (4.2), and the 
other to the horizon region, with coordinates, 



ds 2 = T'dr 2 + A'dr 2 + B'da 2 + 2F'drda + S'dn 2 



(4.7) 



where the charts are related by x = (r + r c ) cos a, y = (r + r c ) sin a and we have taken 
r c = 0.2 for the results shown here. Figure 1 shows the geometry of this chart. Second order 
interpolation is used to translate the metric from one chart into the other. The majority of 
results shown later are for boundary radius R = 3L, which is sufficient to give good results for 
the size of black holes we are interested in. Data is shown mainly for the resolution 80 x 240, 
the first number being the points in the circle direction x, the latter the points in the radial 
direction y. For the localized solution the horizon chart has 80 x 80 points respectively. For 
convergence testing we also consider 20 x 60 and 40 x 120 lattices. The resolution 80 x 240 
is still quite modest, and we expect to improve this in future work. However, it is interesting 
to note that already our results, particularly for the localized black holes, are comparable in 
reach to those performed using other methods at higher resolutions such as in [59]. 

At the asymptotic boundary, chosen at constant y, we have fixed the induced metric by 
specifying the metric functions T, A, S but must still determine boundary conditions for B 
and F. As discussed earlier we impose the geometric boundary conditions £ x = £ y = which 
yield conditions for B and F. Although the charts give the illusion of having boundaries, 
for example at y = and y = R, the metric represented is a smooth Riemannian manifold. 
Hence the boundary conditions of the Einstein-DeTurck equations in a given chart are simply 
determined from requiring this smoothness. On the boundaries x = and x = L one can show 
that all metric functions have Neumann boundary conditions, except F which must vanish. 
At the axis of symmetry of the localized black holes again the functions must be Neumann, 
except F which again vanishes, but in addition regularity of the geometry enforces that S 
vanishes as S = By 2 . On the horizon the metric functions must again be Neumann, except 
for T which vanishes as T = n 2 y 2 A for strings and n 2 r 2 A for localized solutions, where k is 
a constant over the horizon and is related to the temperature as (3 = 2it\/T 00 / n where 
gives the asymptotic value of T. 

There are two classes of behaviour asymptotically in Kaluza-Klein theory. The modes 
with non-trivial harmonic dependence in the circle direction, which behave as massive fields 
from the point of view of the dimensionally reduced theory, decay exponentially at large 
radius. The modes that are constant on the circle, which reduce to light lower dimensional 
degrees of freedom, decay with a power law asymptotically. One subtlety is that taking 
R = 3L for the asymptotic boundary is far enough for all the 'massive' modes to decay away. 
However, for the 'massless' modes which determine the asymptotic value of L and (5 the values 
they reach at R = 3L are not very good approximations to their asymptotic values. Thus 
we find solutions using our domain, applying asymptotic boundary conditions at R = 3L by 
specifying T and A there - let us call these T asym , A asym . We then take the solution near 
the boundary, which is essentially independent of x, and use this as initial data for the ODEs 
governing the massless modes which we then evolve out to a much larger R = 1000L to ensure 
that the massless modes are in their asymptotic regime. Then we extract T^, Aoo from the 
limiting values of T and A there, and in addition the coefficients of the \jr terms in order 
to compute mass and entropy. The difference between , A^ and our control parameters 
T aS ym, A asym is small, of order a few percent, but for accurate determination of physical 
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quantities it is important to perform this 'asymptotic improvement'. Since this difference is 
so small the actual asymptotic values of T, A differ little from their values specified by the 
control parameters - if this were not the case it would be difficult to specify the solution 
of interest by fixing the control parameters. Checking the Smarr law is satisfied then gives 
a gauge as to whether this asymptotic extraction of the charges is functioning correctly. A 
similar procedure was used in [30]. 

Since we are exploiting the various isometries of the metric in order to reduce the effective 
dimension of the problem, the equations contain naively singular terms such as (S — By 2 )/y 2 
at the axis of symmetry y = 0, which derive from the curvature of the sphere which shrinks 
there. Analysis of the equations shows that in addition to the smooth behaviour we require 
there are also singular behaviours. The previous approach to solving the static Einstein 
equations in this context suffered from developing the singular behaviour unless great care 
was taken. However, in both the Ricci-flow and Newton method approach we use here we 
find a very robust behaviour at the axis of symmetry, with no instability issues arising there. 
For both methods, since the approach manifestly is based on the geometry being smooth, the 
potential singular behaviour seems not to surface. For Ricci-flow one can understand this as 
the Ricci-flow attempting to locally diffuse away curvature. 

In Newton's method we must construct the linearization of R H to give the operator Ah 
and solve the linear system A H l R H . It is crucial for the method that the ability to invert 
the operator Ah does not require it to be positive. In practice after the finite differencing 
we map the interior points of our charts into elements of a finite dimensional vector. Then 
we may represent the metric and its curvature R H as vectors. The linearization of R H is 
then computed by a numerical functional differentiation of R H with respect to the metric, 
to yield the matrix operator Ah- Care is taken to implement the boundary conditions and 
interpolation from other patches during this numerical differentiation. Since we have used real 
space finite differencing, the operator Ah is very sparse. Then solving the finite dimensional 
linear system A H l R H can be efficiently achieved using the biconjugate gradient method. In 
practice this stably solves the sparse system, and being based on Krylov methods which utilize 
the finite dimensional nature of the linear problem rather than any form of iteration, is totally 
insensitive to whether eigenfunctions are positive or negative. Only zero modes will present 
problems for this method, and near zero modes may potentially destabilize the method. 



4.4 Ricci-flow applied to localized black holes 

We now demonstrate use of the Ricci-flow method to find small localized 5-d black holes. 
Whilst the Ricci-flow method can be made to work in practice, as we shall see it is less 
elegant than Newton's method if one simply wants to find solutions. The first task when 
constructing these solutions is to choose a suitable one parameter family of smooth initial 
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data that pierces S. We have used the following family, 
ds 2 = (l(r + r c ; d m i n , d max) 9luv 

+ (1 - I(r + r c ; d min , d max )) 92^) dx^dx' v 
g lflu dx ,fM dx' v = T dT 2 + dr 2 + (r + r c ) 2 {da 2 + sin 2 adttl) 
= T dr 2 + dx 2 + dy 2 + y 2 d£l\ 

g 2av dx^dx v = r -^dT 2 + dr 2 + (r 2 + r 2 ) {da 2 + sin 2 adfi|) (4.8) 

which gives a linear interpolation between the near horizon behaviour g 2 and the flat metric 
gi, where we use the interpolation function, 

I{x;x min ,x max ) = -tanhtan ir . (4.9) 

^ \ %max %min J 

We note that this data satisfies all relevant boundary conditions, with k = I/vq. The family 
possesses two parameters, Tb,ro. One combination determines the inverse temperature as 
(3 = (2-7T / ' k)^/Tq. This is then fixed by the boundary conditions, and we are left with one 
parameter to vary. Let us take this to be ro. 

Consider a value of (3 where a small localized solution exists. Since these solutions possess 
one negative mode, the surface £ is codimension one, and we must identify if we are on the 
'+' side or '-' side. It is a question of geometric interest where the '+' and '-' flows from the 
fixed point go. We might expect the '-' flow to behave in a similar manner to that of the usual 
Schwarzschild solution and this is a good guess. Indeed the horizon shrinks to zero size in 
finite flow time, and in a singular manner (with curvature invariants such as the Ricci scalar 
becoming singular there). Presumably if one resolves the singularity the flow continues to 
hot Kaluza-Klein space. More subtle is the what happens to the '+' flow. For Schwarzschild 
in a finite cavity the flow goes to the large black hole. In the asymptotically flat case the 
flow goes in the same direction as the cavity but eats up all of the space, as there is no fixed 
point to reach. However in the case of a localized black hole, in a cavity with boundary at 
R » L one would not expect a large localized solution, as it would be highly deformed 
in the Kaluza-Klein circle direction. Instead the relevant end point is likely to be a large 
uniform string. In order to reach this solution the flow would have to change the topology of 
the manifold from the localized to string topology. The likeliest conjecture then is that the 
'+' mode also flows to a singularity in finite time, although this singularity arises from the 
time circle on the symmetry axis y = shrinking to zero size - if resolved this would then 
flow to the large black string. In the simulations this is indeed what we see. 

We simulated the flows using an explicit implementation of the DeTurck flow, and chose 
the background metric g simply to be the same as the initial metric. We have only explored 
the low resolution 20 x 60 as we are not so much interested in the solutions themselves 
here, but rather the mechanics of finding them. In figure 2 we show the horizon area and 
time circle 2-sphere area for 3 flows with (3 = 1.5. We find that tuning ro a critical value, 
r criticai — 0.598123, is required to approach the fixed point. This critical flow is shown, in 
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Figure 2: Behaviour of the area of the horizon (left) and area of the time circle 2-sphere (right) as a 
function of Ricci-fiow time for 3 initial data. In blue, r is slightly greater than the critical value, in 
red, ro is slightly less, and the yellow gives the critical value. The red and blue curves diverge from 
the yellow due to the presence of the single negative mode. The dashes horizontal lines give the value 
of the area at the fixed point, calculated using Newton's method, and thus we see agreement with the 
critical flow. 

addition to two other flows, one with a lower tq = 0.598 where the horizon sphere shrinks to 
zero size in finite time - the '-' behaviour - and the other for ro = 0.600 where the time sphere 
shrinks and the horizon expands - the '+' behaviour . We see the flow ro = r cr iti ca i appears to 
stabilize at a fixed point. Interestingly for the '-' and '+' behaviours when the sphere shrinks, 
the time circle expands, and vice versa. In figure 3 we show full proper embeddings of the 
horizon and symmetry axis against flow time for these three flows, and we see that the entire 
horizon remains round as it shrinks while the symmetry axis length seems to be unaffected. 

Since the nature of the two possible singularities is straightforward to determine, it is 
simple to automate a bracketing algorithm that finds the solutions. In figure 4 we show 
brackets in ro that include the value of r cr i t i ca i against (3 for a range of (3. The value of 
ro in the initial data gives the initial size of the horizon. In the plot we also show the 
size of the horizon of the fixed point. We see that as (3 increases, the size of the horizon 
increases relatively more quickly than the size of the horizon of the fixed point. Indeed we 
see that for quite small (3 ~ 1.6 the initial horizon diameter (ie, 2ro) is already becoming 
comparable to the size of the Kaluza-Klein circle. We find that to accurately simulate the 
flow we require higher resolutions, and it appears that the initial data used may develop a 
coordinate pathology during the flow. Presumably this could be resolved with more effort, 
but we have no attempted to do so here. Even doing this, it is not clear how high (3 could 
be pushed while allowing a critical solution - recall that our initial data does not have to 
intersect T, for all 0. We see in the next section that Newton's method also doesn't find 
localized solutions above (3 ~ 3.39. Thus the range of (3 we have probed here is substantially 
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Figure 3: Cross-sections of proper embeddings of the horizon are shown in red for L — 1 against 
Ricci-fiow time for the same 3 flows as in the previous figures. The axis of symmetry is also shown 
in blue with its proper length. The top figure corresponds to the red curves, r less than the critical 
value, the middle figure gives the critical flow, and the bottom figure corresponds the blue curves so 
ro is greater than the critical value. 

less than that which is possible, and given the growing size of ro relative to the fixed point 
size, we suspect that our family does fail to intersect S before this maximum (3. 

We have only attempted to find the small localized solutions with their single negative 
mode. Whilst in principle we believe that one can find the large localized solutions and non- 
uniform strings with two negative modes by finding an appropriate 2 parameter family of 
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Figure 4: Brackets for r$ over a range of f3. For the lower end of the bracket the sphere collapses, 
and for the upper end the time circle shrinks in finite flow time. Outside this range, higher resolution 
is required to accurately simulate the flows. The value vq gives the initial radius of the horizon for 
the flow. We also plot the equatorial radius of the horizon for solutions shown as small dots, together 
with the small localized approximation for this radius which is rather accurate over this range of [3. 



initial data, we have not tried it. The principle complication is identifying how to tune the 
initial data, given the late time behaviour of the flow, and it would be interesting to examine 
whether there is interesting geometry involved, as in the one negative mode case, or whether 
this becomes purely a technical problem. 

Needless to say the fixed points found agree with those found by Newton's method - 
for example the areas shown for the critical flows in figure 2 agree with the values found 
from Newton's method - and hence we defer showing the properties of the solutions until the 
next section. We see the Ricci-flow algorithm we present does work in practice. It is simple 
to implement the explicit flow, and it would be relatively straightforward to implement an 
implicit algorithm to allow larger time steps to be taken than in the explicit case although 
whether this would improve the speed is unclear as this depends on the details of solving the 
implicit system. 

Whilst the flow itself is easy to construct, and the tuning process easy to implement, at 
least for codimension one E, it appears that the tuning must be done rather precisely to closely 
approach the fixed point. This is to be expected due to the asymptotically flat boundary 
condition. The approach to a stable fixed point is governed by the lowest eigenmode. We 
truncate our asymptotic region to have a boundary at finite radius. Then we expect the lowest 
mode is goverened by that radius. Taking the radius to be larger will reduce the eigenvalue 
of the lowest mode, and removing the truncation will yield a continuous spectrum down to 
zero eigenvalue. In practice the size of the bracket must be systematically reduced until the 
flow approaches and remains near the fixed point for some considerable flow time. In our 
case, for R = 3L, tuning ro to 10~ 6 yields flows that remain near the fixed point only for 
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1.0 < A < 2.0 before it develops a '-' or '+' singularity. Taking the solution around A = 1.5 
then gives a good approximation to the fixed point solution, but this is an additional source 
of systematic error beyond the finite differencing and infinity truncation, and removing the 
infinity regulator makes this problem worse. 

Thus whilst conceptually and practically simple to implement, the Ricci-flow method does 
suffer from the problem that it is difficult to attain good accuracy of the fixed point solution 
by tuning a family of initial data. For low resolutions such as the 20 x 60 used here the error 
in the resolution and finite infinity truncation may outweigh this accuracy error. However, for 
higher resolutions, and for removing the regulator clearly this is presents a serious deficiency 
of the algorithm that limits its usefulness. 

Such a precision tuning problem could be avoided by modifying the algorithm as follows. 
A crude bracket is constructed using the initial family of data above. Each side of the 
bracket flows close to the fixed point for some time and then diverges away in the direction 
of the negative mode. However, we may take the metrics given by flowing each side of the 
bracket for a time so that the flow has approached the fixed point, and not yet diverged away. 
Let us denote these gi^- Then one may construct a new one parameter family of metrics 
g{\) = A <7i + (1 — \)g2 which again bracket the fixed point for A € (0, 1). One then proceeds 
to use this new family to refine the bracket. Since one is closer to the fixed point, for a given 
precision of tuning in the parameter A one should get closer to the fixed point than for an 
equivalent tuning of the initial family parameter r$. Clearly this procedure can be iterated. 
Whilst this procedure will work, it is somewhat ad hoc. Deciding how long to run the flow 
for, and how accurately to find the bracket, are just some of the points that one would have 
to specify. 

In essence the Ricci-flow can get close to a solution by the simple tuning of an initial 
family of data, but once one is near the fixed point, there are more efficient methods to hone 
closely in on that fixed point. Newton's method described next represents precisely such a 
method of honing in on the fixed point, and thus one use of the Ricci-flow method could be 
to provide initial data for Newton's method, where this was difficult to construct. In the case 
of the solutions we discuss here we have analytic initial data in the basin of attraction of the 
solutions of interest for Newton's method, and once one solution on a branch is found, one can 
easily move along the branch. However in more complicated situations it might be unclear 
how to write down initial data in the basin of attraction of a branch of solutions, particularly 
if this branch were disconnected from other known branches. Then Ricci-flow might provide 
a valuable tool to construct metrics close to the fixed point. 

We note that the Ricci-flow, being a flow in the space of diffeomorphism classes, does 
give global geometric information about the space of metrics as the flow proceeds. Newton's 
method, while simpler to use in practice, gives no geometric information other than the 
solution it may find. For example, suppose we were not interested in the details of a solution, 
but merely its existence - black holes on Randall-Sundrum branes are an interesting problem 
of this variety [98-100]. Then a relatively coarse bracketing of the '+' and '-' behaviours gives 
strong evidence of existence, even though the flows would not go very near the fixed point 
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and therefore the details of the fixed point solution could not be read off. We may regard 
the results above as giving a 'numerical proof of existence of localized black hole solutions 
- although in this case there is no reason to doubt their existence. With Newton's method 
one must test existence by studying convergence with discretization scale. However with the 
Ricci-flow, since we are clearly probing the global structure of the space of geometries, rather 
than just finding a single solution, we may regard our brackets above as providing convincing 
evidence for existence. 

As we will see shorty, in the case of the localized solutions Newton's method clearly 
demonstrates that the small localized branch reach a maximum (3 which can be continued 
to large localized branch, so whilst (3 has a maximum the solution branch is smooth there. 
The Ricci-flow method is compatible with this in the sense that we have not found solutions 
for higher (3. However, one might imagine a branch of solutions could actually end at some 
(3. In this case it might be difficult to know from Newton's method why solutions were not 
being found. With the Ricci-flow method, one might be able to gain some additional global 
information from the nature of the flows. 

4.5 Newton method applied to localized and non-uniform string solutions 

The implementation of Newton's method is technically considerably more complicated than 
the Ricci-flow, as the Einstein-DeTurck equations and metric must be recast as vectors, the 
operator Ah must be computed - we use a numerical functional differentiation - and then 
the large but sparse linear system must be solved - in our case using biconjugate gradient - 
in order to provide the basis for Newton's method. However, once this technology is in place 
the method works extremely well in practice. 

The first step in finding the solutions of interest is to find appropriate initial data in 
the basin of attraction of that solution. This is not trivial, but in practice is manageable 
with a little trial and error. For example, in the case of 5-d localized solutions, we find that 
the initial data given by equation (4.8) above with r$ = 0.15 gives a good initial metric and 
background metric to find the (3 = 1.2 small localized solution. Once a solution on small is 
found it is very simple to find neighbouring ones simply by perturbing the solution metric and 
background metric to change (3 to a new nearby value. For example, we use the perturbation, 
5g = I(y, L,2L)dr 2 to implement this. For a small perturbation we remain in the basin of 
attraction of the solution at the new value of (3, and that solution is quickly found. 

A subtlety then arises when a maximum/minimum in j3 arises along a branch of solutions. 
For the example of small localized 5-d black holes we start with (3 = 1.2 and then gradually 
increase (3 but find no convergence above (3 = 3.39. However, plotting physical quantities 
against (3 one sees that this is not the end of the branch of solutions, but merely a maximum 
in (3 so that (3 max = 3.39. In order to pass this maximum to find the large localized solutions 
one must again find initial data with (3 < f3 max but that is in the basin of attraction of the large 
localized branch rather than the small localized one. This may be achieved by taking solutions 
found for small branch near the suspected maximum (3 max , and using them to extrapolate to 
an initial guess with (3 less than the maximum which should lie on large branch. Suppose we 
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have small branch solutions gi, with i = 1, ... N, and corresponding such that < /3j if 
i < j. In order to extrapolate we must find a quantity which varies monotonically through 
ftmax, rather than having a maximum/minimum there. For example, we may take the quantity 
Laxis which appears to decrease monotonically along the small branch near (3 max and across 
the join to the large branch. Thus at the maximum, (3 varies quadratically with L ax i s . Then 
we use the values of (L ax i s )i for the solutions, and extrapolate the data to find an initial 
guess gN+i for a value (L ax i s )N + i which is small enough to correspond to a large localized 
solution. (3n+i is then simply determined from gw+i- In practice this procedure works well 
and gives good results, allowing us to cross to the large localized solutions, having found the 
small localized branch. 

Finding initial data for the non-uniform strings is less straightforward than for the small 
localized solutions. It is very straightforward to find uniform strings, using an analogous 
scheme to the small localized solutions - simple analytic initial data can be found in the 
basin of attraction. In particular in 5-d we use, 

2 2 

ds 2 = = -^-^dT 2 + dx 2 + dy 2 + (r 4 (l - I(y, 0,1)) + y 4 ) 1/2 d£l 2 (4.10) 
y 2 + r 2 

to find uniform strings, and choosing ro = 0.28 yields a solution near the marginal point 
which is at [3gl = 1-752. The challenge however is to find initial data which is in the basin 
of attraction of the non-uniform strings. In 5-d we know from Gubser's work [87] that for 
weakly non-uniform strings we have (3 > (3gl- Thus we perturb the metric and background 
metric for this marginal solution to lower its temperature, 

SA = A(l - 0.02(1 - I(y, 0, L))) (4.11) 

and perturb only the metric (not the background) with a non-uniform profile, 

5B = B (1 - 0.05 cos ttx (1 - I(y, 0, L))) 

SS = S (1 - 0.05 cos ttx (1 - I(y, 0, L))) (4.12) 

We find this does indeed relax to a weakly non-uniform string. We may then traverse the non- 
uniform branch by adjusting the temperature as for the localized black holes. Note that it is 
convenient to keep the background metric uniform as then it is straightforward to distinguish 
between the uniform and non-uniform branches by simply looking at whether the metric is 
uniform. 

However a more systematic method is to find the uniform string at (3gl, an d construct 
the Gregory-Laflamme zero mode. This is simple in practice as Newton's method already 
requires the linearization A# to be computed, and it is straightforward to compute the low 
lying eigenfunctions of a sparse matrix such as this. Having done so, we take a uniform 
solution with (3 > (3gl but still near the marginal point, and then add some amount of the 
zero mode to it. We note that the linear perturbation of the (3 = (3gl, marginal solution 
preserves most of the boundary conditions even when added to a solution with a different 
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P > Pgl- 8 Adding very little, one remains in the basin of attraction of the uniform solution. 
However, when sufficient amplitude of this perturbation is added, one is moved into the basin 
of attraction of the non-uniform solution. Having found one non-uniform string, one may 
explore the branch by perturbing the temperature as for the localized solutions. 

A potentially attractive feature of Newton's method is that in principle the convergence is 
uniform in the wavelength of the deformation from the fixed point. Removing the asymptotic 
truncation for Ricci-flow presents the problem that the long wavelength modes evolve slowly, 
whereas for Newton's method this is not the case. However in practice the time does depend 
on how long the bi conjugate gradient method takes to find solutions, and this can depend 
on the spectrum of A#. We have not investigated this here, and other methods to solve the 
linear system may behave differently, but it is a point to note. 

In figures 5, 6 and 7, we show various how physical quantities vary against (3 for the 5-d 
localized and non-uniform strings. Data is shown for the highest resolution used, 80 x 240. 
Where available we give the perturbative results for localized solutions (computed up second 
order [82]) which compare well up to (3 ~ 2 where the perturbative expansion at different 
orders starts to differ significantly, signaling the non-perturbative regime. The limited data 
for 5-d localized solutions previously existing agrees [32,59], and the previous calculations of 
5-d non-uniform strings also agree [39,40]. In appendix A we estimate the numerical error in 
the results to be of order ~ 2% or less by considering data computed at different resolutions 
and for different boundary positions. In that appendix we also check that the solutions are 
consistent with being Ricci-flat rather than non-trivial Ricci solitons. 

It was previously observed for 6-d that such physical quantities were consistent with a 
merger of the localized and non-uniform branches [59]. Our 5-dimensional data is indeed 
consistent with such a merger at j3 m erger — 2.64, where the behaviour of various quantities 
plotted is consistent with being continuous, although not smooth. The figures show agreement 
in all quantities, taking into account the percent level systematic errors. We find experimen- 
tally that increasing the resolution allows solutions to be constructed closer to this potential 
merger point, with 80 x 240 being sufficient to approach quite closely. We note that we are 
able to find localized solutions closer to the merger point than previous numerical approaches 
which employed higher resolutions [59]. 

In figure 5 we show the various quantities related to thermodynamic behaviour. The free 
energy curve shows the expected first order phase transition between small localized solutions 
at high temperature, and uniform strings at low temperature, with the transition being at 
(3 ~ 2.72. The free energy shows a clear cusp behaviour at P max - The mass, area and tension 
appear to have a smooth behaviour at f3 m ax and we return to this interesting point later. In 
figure 6 we plot various geometric properties of the solutions and in figure 7 we show geometric 
quantities that characterize distance to a possible merger. For the non-uniform solutions 
distance from the merger can be characterized by the minimum radius of the horizon, R m in- 
The smallest value this reaches for our data corresponds to A ~ 3.6, where A is the quantity 

8 Only the boundary condition for £ = at the asymptotic boundary is not preserved, and this is quickly 
restored during Newton's method. 
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Figure 5: Thermodynamic behaviour of the various static Kaluza-Klein black holes for L = 1. Top 
left: Plot of mass, M, against j3. The small localized branch are shown in red, the large localized 
branch in orange, and the non-uniform solutions are in light blue. All are for 80 x 240 resolution 
with R — 3L. The uniform branch is given by the dark blue line, and the small localized analytic 
approximations are shown by a thin line. The solid blue disc indicates the Gregory-Laflamme marginal 
point. Similar conventions are used to plot the horizon area A^oriz (top right), tension r (bottom left) 
and free energy F (bottom right). 

used by Gubser to characterize the non-uniformness, A = (R max /Rmin — l)/2 [87]. On the 
localized side we may assess the distance from merger either by measuring the proper length 
of the symmetry axis L ax i s , or by measuring the maximum radius of the Euclidean time circle 
on that axis which, for these solutions, is given by R T . As we see from the figure both these 
shrink, where we have plotted R ax is = L ax i s /2ir so that this may be treated as a radius and 
all three quantities - R m in, Ra X is and R T - may be directly compared in value. We see that 
the localized branch is not found quite as near to the merger point as the strings, but still 
they are close. We note that R ax i s ~ Rt as they decrease, indicating that the 2-sphere formed 
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Figure 6: Aspects of the geometry of the various static Kaluza- Klein black holes. The same conven- 
tions are used as in previous figure. Top left: Plot of equatorial horizon radius, R eq for the small and 
large localized solutions and maximum horizon radius R m ax for the uniform and non-uniform strings 
against f3. Top right: Plot of polar horizon length, L po i ar , for the localized solutions and horizon 
length, Lhoriz, for the strings against (3. Bottom: Plot of area eccentricity, e, for the small and large 
localized solutions against (3. 



by the fibration of the time circle over the symmetry axis is becoming round as it shrinks, 
which we would expect in the type of merger discussed by Kol [88] . 

In figure 8 we show a plot of the horizon embeddings as a function of (3, and again see 
consistency with a merger, where we recall that there is percent level systematic error in the 
geometry of the embedding, and the value of (3. It would be interesting to attempt to check if 
the conical geometry proposed by Kol [88] near the merger point is emerging on the localized 
side, as has been checked on the non-uniform side [39,93]. We leave this for future work. 

Since we are constructing the linearized operator Ah, it is straightforward to find the 
low lying eigenvalues and functions of this sparse matrix. However, in practice we have 
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Figure 7: For the small and large localized solutions we show the axis length, L ax i S , and mirror 
plane time circle size, R T , against (3. For the strings we plot minimum horizon radius R m i n - These 
curves characterize how near to the suspected 'merger point' the solutions constructed reach. The 
same conventions are used as for previous figures. 

utilized the isometries of the solution, and so only compute the operator Ah restricted to 
action on modes that are constant on the time circle and 2-sphere. Hence we may only find 
eigenmodes that are singlets under the U{\) x SO(3) isometry group. In principle a good 
way to characterize the merger is to compute the spectrum of this operator and examine 
its behaviour on both sides of the merger point, as the low lying spectra provide a probe of 
the geometry globally. Firstly we compute the negative eigenvalues showing these in figure 
10. As noted earlier in section 3.1 these negative eigenvalues should correspond to negative 
eigenmodes of the operator A^. We see that as expected the small localized solutions have 
one negative mode in the class of U{1) x SO(3) invariant eigenfunctions, and we find that the 
large localized solutions have two. At the join of these branches, the tangent to the space of 
solutions will precisely give a normalizable (ie. in our case, boundary condition preserving) 
zero mode, and it is this which continues to the second negative mode for the large localized 
branch. We find the non-uniform strings have two U(l) x SO(3) invariant negative modes. In 
principle additional negative modes might develop nearer to the merger point than we have 
accessed. However it is interesting that for the range of solutions explored here we see two 
U{\) x SO{3) invariant negative modes as the merger point is approached from both sides. 
In particular a new result is that over the range of (3 tested we see no zero modes develop 
on either branch that would correspond to that branch joining with some exotic new branch 
of static solutions with U(l) x 50(3) isometry. The only zero modes are the one at the 
Gregory-Laflamme marginal point, and that at (5 = (5 max which is understood simply as a 
maximum in (3, and not as a mode that may be exponentiated to generate an independent 
branch of solutions. We see that one negative mode appears to become increasingly negative 
near the merger point and appears to be diverging to minus infinity - presumably there is an 
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Figure 8: Cross-sections of proper embeddings for the horizons of the various black holes against j3 
for L = 1. The upper frame shows the small localized solutions, and the lower frame the large localized 
and non-uniform string solutions. The horizon is shown in red, and for the localized solutions the axis 
of symmetry is also shown with its proper length in blue. 



associated critical behaviour that might be studied. The other negative eigenvalue appears 
to remain finite, and is consistent with a merger. 

In figure 11 we show the lowest 30 positive U (1) x SO(3) invariant eigenmodes for the large 
localized and non-uniform solutions near the merger for R = 3L. These modes correspond 
to physical eigenmodes (as opposed to short wavelength modes which are sensitive to the 
discretization) of the various operators discussed in section 3.1, namely Ah and Ay- Without 
the asymptotic boundary regulator this spectrum would be continuous down to zero. However 
with the regulator R = 3L it becomes discrete and can be used as a global probe of the 
geometry, and in particular can be used to test for the potential merger. A minor point is 
that due to the finite boundary truncation and discretization a few of the modes have a small 
imaginary component in their eigenvalue, and thus it is the real part of the eigenvalue that 
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Figure 9: This figure shows another view of the previous embedding figure. 



we are plotting. We estimate the error in the eigenvalues in this figure due to discretization to 
be quite small, of order ±1 in absolute value, from comparison between different resolutions. 
We see a good correspondence of the spectra at the potential merger point, which again lends 
considerable weight to the idea that these branches merge. We note that assuming a merger 
occurs, it appears unlikely that any of these positive eigenmodes have the opportunity to 
become negative closer to the merger point than we have probed on either side, although of 
course this possibility cannot be ruled out. 

In order for the small and large localized branches to join smoothly - and there is no 
reason to expect this join not to be smooth in the space of metrics - there must be a maximum 
in the mass at some j3 < f3 max . The question is whether the maximum lies exactly at (3 max 
or below it. This is a very interesting question as if the maximum lies exactly at (3 max then 
the specific heat capacity of all the solutions shown will be negative. Note that this would 
mean the mass against /3 relation would be cuspy at j3 max - as for example for the free energy 
- but this does not imply a lack of smoothness in the space of metrics. However, if it lies 
below it, then for some range of (3 the localized solutions may have positive specific heat, as 
the gradient of energy with temperature - ie. mass with inverse (5 - will be positive. From 
the previous figures we see that it is a close call. In figure 12 we reproduce the area, mass 
and free energy plots focussed near the join of the small and large branches. The mass and 
area appear to be smooth at f3 max , whereas the free energy, shown for contrast, clearly shows 
a cuspy behaviour. 

From a purely geometric point of view it would have been surprising if the maximum in 
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Figure 10: The negative eigenvalues, ui 2 , of the operator Ai which are U(l) x SO(3) invariant are 
shown plotted against p for small and large localized, and non-uniform string solutions. The solid 
discs show the results for a uniform string at the marginal Gregory-Laflamme point, where we have 
an zero mode and one negative mode. We observe that the zero mode continues to a negative mode. 
The line gives the leading order approximation for small localized solutions derived from the negative 
mode of 5-d Schwarzschild (computed in [101]). 

mass coincided exactly with a maximum in f3. One might wonder then why this is true for 
free energy? This follows from the fact that the free energy is related to the action. At the 
maximum (5 the tangent to the space of solutions precisely gives a normalizable zero mode. 
Note that the tangent to a solution branch always gives a zero mode, but not a normalizable 
(or boundary condition preserving) one, as it will generically result in a perturbation of 
[3. Since the Ricci-flatness condition derives from the variation of the Euclidean Einstein- 
Hilbert action (plus boundary terms), /, this implies that I should be stationary at (5 max as 
moving along a solution branch will change / only by terms involving the change in boundary 
conditions. It is precisely a normalizable zero mode that leaves I invariant. Since the free 
energy F is simply related to I by I = (3F for these solutions, we should expect it to be 
stationary at (3 ma x ~ resulting in a cusp behaviour for F against (3 - and this is indeed what 
we see. However there is no analogous argument that the mass should be stationary. 

We see that actually the maximum mass lies on the small localized branch a little below 
the maximum j3 point. Let us denote this j3 = I3 mass . Then we find that (3 mass / ' (3 max ~ 
0.995±0.001, where the error is estimated by comparing data for 20 x 60, 40 x 120 and 80 x 240 
resolutions, each of which sees the same effect, and from comparison of R = 3L, R = 6L and 
R = 10L data, which again all see this effect. This maximum mass point coincides with the 
maximum area point, as dictated by the First Law. Whilst the effect we see is robust to 
varying resolution and boundary truncation, the (3 range from j3 mass to f3 max is rather small, 
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Figure 11: The 30 lowest positive U(l) x 50(3) invariant eigenvalues of Ah for the large localized 
(left of the axis) and non-uniform solutions (right of the axis) in the region of the potential merger 
for asymptotic boundary truncation R = 3L. This truncation discretizes the spectrum, which can 
then be used as a probe of the potential merger. For both localized and non-uniform solutions the 
eigenvalues are coloured in the order given by their value. We observe that there appears to be a 
good correspondence between the spectra of the localized and non-uniform solutions near the merger, 
consistent with the spectra being continuous (although presumably not smooth) across the merger. 

and so one might worry whether our numerics are correctly resolving this detail sufficiently 
well. That said, it does appear that the curves for mass and area against (3 - in contrast to 
the free energy - look rather smooth, and this alone implies (3 m ass < Pmax- Clearly this is an 
issue for future work to confirm. However, if confirmed, this narrow window of small localized 
solutions, from f3 ma ss to Pmax would represent the first examples of static black holes that have 
positive specific heat, but contain normalizable negative modes. Whilst it has been shown 
that there is a relation between local thermodynamic instability (for static vacuum solutions, 
a negative specific heat) and the existence of negative modes [53-56], the solutions exhibited 
here represent a counter-example to any argument that suggests the converse might be true. 
This would be an interesting complement to the recent analytic arguments that AdS-Kerr 
provides a stationary counter-example [62] . 

5. Summary 

We have provided a new general approach to solving static and Euclidean Einstein equations 
in vacuum. Whilst we haven't included matter, the methods should generalize trivially, as the 
stress tensor terms will not change the ellipticity property of the Einstein-DeTurck equations, 
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Figure 12: Plots of horizon area (top left), mass (top right) and free energy (bottom), which are 
reproductions of previous figures but focussing on the join between the small and large localized 
solutions. In particular the mass and area curves appear smooth at the maximum (3 point, and 
therefore indicate a positive specific heat for a narrow range of the small localized solutions. In 
contrast the free energy curve is cuspy at the maximum (3. 



and the matter equations can be phrased as elliptic equations too 9 . Then the Ricci-flow 
method naturally generalizes to a diffusion flow of the Einstein-DeTurck-matter system, and 
Newton's method extends trivially. 

Our method works for general cohomogenity metrics, although it allows isometries to be 
used to reduce the effective dimension of the problem. For a stable fixed point our Ricci-flow 
method simply reduces to DeTurck flow, used successfully previously to study Calabi-Yau and 
del Pezzo geometries. It would be interesting to return to these problems and try applying 
our Newton method to compare with Ricci-flow. However, the new feature of our approach is 



Scalar field equations will already be elliptic, and gauge fields can be made elliptic using a similar trick to 
that of DeTurck. 
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its ability to deal with Euclidean solutions with negative modes - either by tuning data and 
using Ricci-flow, or by using Newton's method. We have seen that in the case of Kaluza-Klein 
black holes the latter approach provides the more practical route, and yields excellent results. 
However the Ricci-flow method does have the virtue of being a flow in the space of geometries, 
and can potentially be used to give global information about the space of solutions, and in 
particular a 'numerical proof of existence. 

Our example application to Kaluza-Klein black holes has provided some interesting new 
results in itself; 

• We have seen that the localized solutions have a minimum temperature and hence divide 
into the two branches we term 'small' and 'large'. 

• We have seen consistency with a merger between the localized and non-uniform strings, 
having constructed localized solutions closer to the potential merger point than previ- 
ously found. 

• We have computed the Euclidean negative modes of these solutions which are constant 
on the time circle and 2-sphere, finding the small localized solutions have one, and the 
large localized and non-uniform solutions constructed have two. One of these two has 
finite eigenvalue at the potential merger, the other apparently diverges to minus infinity. 

• Consideration of the low lying positive spectrum of Ah which is constant on the time 
circle and 2-sphere also shows consistency with a merger, and indicates that there are 
no additional zero modes likely near the merger point which preserve the U(l) x SO (3) 
isometry. The existence or not of additional zero modes preserving the isometries is 
important when considering the various merger scenarios [88] and global structure of 
the phase diagram - for example when including multi-black hole configurations [102]. 

• A more speculative but rather interesting result is our evidence for an intriguingly nar- 
row window of small localized solutions with positive specific heat but which possess 
two normalizable Euclidean negative modes. Whilst local thermodynamic instability 
appears to imply the existence of negative modes [53-57], such solutions would then 
represent counterexamples to any reverse claim. We caution that this window is sur- 
prisingly narrow and thus this is a somewhat subtle numerical issue that will require 
future work to confirm with confidence. 

We began by stating that the problem of directly finding general stationary, static and 
Euclidean solutions in a unified framework was a problem of interest. We have managed 
to find a way to unify the static and Euclidean approach for geometries with and without 
a horizon. However our method clearly does not easily generalize to the stationary case, 
where an analytic continuation to a (real) Euclidean geometry cannot be performed. Thus 
an interesting question for future work is how to unify the stationary case with the other two. 
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A. Convergence and boundary tests 
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Figure 13: Horizon area (left) and mass (right) for 3 resolutions 20 x 60 (red), 40 x 120 (blue) and 
80 x 240 (purple), all with R — 3L. Correct second order scaling to the continuum is observed, and 
the 80 x 240 data shown elsewhere in the paper is estimated to have a percent level finite resolution 
error. 

Here we estimate finite resolution error and boundary truncation dependence. In figure 
13 we show the horizon area and mass curves giving data for 3 resolutions; 20 x 60, 40 x 120 
and 80 x 240, all for R = 3L. We see correct scaling to the continuum in these data - and 
all other data shown in the paper - and we estimate the error in the 80 x 240 data relative 
to the extrapolated continuum at approximately the percent level for the various quantities 
shown in this paper. We find that we may obtain non-uniform solutions and large localized 
solutions closer to the possible merger point as we go to higher resolution, and therefore it is 
likely that at a given resolution this error increases near to the merger point. 

In figure 14 we show the effect of changing the finite boundary truncation, giving the 
horizon area for R = 3L (the value used elsewhere in the paper), R = 6L and R = 10L. 
We see that the systematic error in the R = 3L results are less than 2%. We find the same 
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Figure 14: Horizon area for 3 boundary truncations - R = 3L (red), R = 6L (blue) and R = 10L 
(purple), all with resolution 20 x 60. We observe that already for R—'iL the error introduced by the 
boundary truncation appears to be small, of order < 2%. This is true for other observables studied. 
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Figure 15: Fractional violation of the Smarr relation for the localized solutions (left) and non-uniform 
strings (right) for 3 resolutions 20 x 60 (red), 40 x 120 (blue) and 80 x 240 (purple), all with R = 3L. 



behaviour for other measured quantities. The Smarr relation in equation (4.6) allows another 
global check of error. In figure 15 we plot the fractional error in this relation for the localized 
black holes, 



Frac Error 



3TS + tL 
2M 



(A.l) 



and see that for the 80 x 240 resolution this fractional error is less than 1% for all solutions 
found. We note that there are two contributions to error; the discretization error, and the 
finite boundary error. We expect that the discretization error will dominate, particularly 
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for small localized solutions, and the data is consistent with this, with increasing resolution 
decreasing the error. For the 80 x 240 data, we expect that the major source of error is the 
finite boundary truncation, and for the non-uniform strings the error is at a similar level to 
that estimated from the finite truncation above. For the localized solutions the error is less 
than this, but we suspect this is coincidental. 



Figure 16: Maximum absolute value of the vector £ for the small localized solutions for 3 resolutions 
20 x 60 (red), 40 x 120 (blue) and 80 x 240 (purple), all with R = 3L. The left frame shows max \£ x , f v | 
taken over the Cartesian chart, the right frame shows max|£ r ,£ Q | taken over the near horizon chart. 
We are interested in the behaviour of these maximum absolute values as the resolution increases - the 
values themselves are not coordinate invariant quantities and therefore carry no significance. We see 
that at a given /3 the maximum value of decreases in accordance with second order scaling to zero 
in the continuum. 

Finally, in figure 16 we show that the vector field £ is indeed consistent with vanishing 
in the continuum, confirming that we are finding a Ricci-flat solution, rather than a Ricci 
soliton. The maximum absolute value of the components of the vector field £ in the two 
charts is plotted for the small localized solutions for the 3 resolutions 20 x 60, 40 x 120 and 
80 x 240. For all (5 the maximum absolute value is seen to decrease consistent with second 
order scaling to zero in the continuum - the values themselves, being coordinate dependent, 
are not physical. We see the same behaviour for the large localized and non-uniform strings. 
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